Left-Invariant Diffusion on the Motion Group in terms 
of the Irreducible Representations of SO (3) 

Marco Reisert and Henrik Skibbe 

University Medical Center Freiburg, Medical Physics, Germany 
marco. reisert® uniklinik-freiburg. de 

February 27, 2012 

Abstract 

In this work we study the formulation of convection/diffusion equations on the 3D mo- 
tion group SE(3) in terms of the irreducible representations of 50(3). Therefore, the left- 
invariant vector- fields on SE(3) are expressed as linear operators, that are differential forms in 
the translation coordinate and algebraic in the rotation. In the context of 3D image processing 
this approach avoids the explicit discretization of SO (3) or S2 , respectively. This is particular 
important for SO(3), where a direct discretization is infeasible due to the enormous memory 
consumption. We show two applications of the framework: one in the context of diffusion- 
weighted magnetic resonance imaging and one in the context of object detection. 

Keywords: Spherical Harmonics, Left-invariant Diffusion, Partial Differential Equations, Wigner D-Matrix, Clebsch- 
Gordan coefficients, Diffusion weighted MR-imaging, High Angular Resolution Diffusion Imaging (HARDI), Diffu- 
sion Tensor Imaging (DTI), Spatial Regularization, Spherical Hough Transform 

1 Introduction 

Image Processing in 3D becomes more and more popular and necessary due to the enormous 
amount of scientific data acquired with modern imaging techniques like magnetic resonance imag- 
ing, computer tomography or confocal laser microscopy to mention only a few. Still, the processing 
of directional or tensorial information derived from primary modalities or directly measured like 
in diffusion weighted imaging (DWI) drives todays hardware to its limits. For example, imagine a 
typical DWI measurement with imaging matrix of 100 3 . If we want to sufficiently represent the ori- 
entation space with e.g. 500 points, we need already about 4 GB of memory for just one instance. 
As typical algorithms (e.g. like Conjugate Gradients) usually require for more than one instance, 
we are already at the limits of a common desktop PC. In this article we describe how the generators 
of diffusion and convection on R 3 x S 2 and R 3 x SO (3) can be described and implemented in terms 
of the irreducible representations of the 3D rotation group. Most of the implementations solving 
R 3 x £ 2 -diffusion equations E El El HI relied on an equiareal discretization of the two-sphere 5 2 , 
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while implementation for functions SE(3) \-> C even do not exist due to the enormous memory 
consumption, although there might be several useful applications like feature detection for non- 
rotation symmetric templates. Indeed, there are implementations [3 ] that use spherical harmonics 
as an intermediate ^-interpolation scheme, but they cannot benefit of the well-known advantages 
of the spherical harmonic representation, like the compact and memory efficient storage, the ana- 
lytic and efficient computations of ^-convolutions, and the closeness under rotations. The aim of 
this paper is the formulation of common differential operators acting on functions M 3 x S 2 H> C 
in terms of spherical harmonics. More generally, we show how a large class partial differential 
equations on functions SE{3) \-> C can be represented in terms of the irreducible representations 
of the rotations group SO (3) and solved without any angular discretization. In this way one can 
benefit from all the advantages which the harmonic representations offer. A discretization of the S 2 
or SO (3), respectively, is avoided, and one is able to implement diffusion on the full SE(3) with 
reasonable memory consumption. We show one application in the context of high angular resolu- 
tion diffusion MR-imaging, where spherical harmonic representations are common and compare 
to the equiareal representation. Second, it is shown how the proposed framework can be used to 
implement the spherical Hough transform in an efficient manner. 

1.1 Related Work 

In the context of line and contour enhancement in 2D the special motion group SE{2) plays a 
key role [[H [71 [8). It can be used to set up a scale space theory. More recently, extensions 
to 3D of these concepts appeared [0. While the applications in 2D are typically related to feature 
detection and image enhancement, the 3D extension offers a new application field: the processing of 
diffusion weighted magnetic resonance images (DWI). In DWI already the acquired measurements 
are functions on K 3 x 5 2 . Based on the directional dependency of water diffusivity in fibrous 
tissue of the human brain it is possible to reveal underlying connectivity information. One of 
the main challenges in DWI is the estimation in so-called fiber/diffusion orientation distributions. 
There are numerous methods for estimating orientation distributions: classical Q-ball imaging [9 J, 
constrained spherical deconvolution [10], proper probability density estimation (TTl[T2l[T3l[T4| and 
spatially regularized density estimations for tensor- valued images Il5l|4l[T6l[T7l[l8l[l9l. Most of 
the employed algorithms rely on tensorial or spherical harmonic representation of the orientation 
distributions. On the other hand, most of the algorithms for orientation distribution estimation that 
consider the local surrounding of a voxel and use intervoxel information rely on a discretization of 
the two-sphere flU El El SJ . 

In two dimensions the representation of orientation and tensor fields in terms of circular har- 
monics (or, the irreducible representations of 50(2)) is relatively simple and quite frequent in liter- 
ature (20l|2Tl[22l[23]|. Complex Calculus offers a well-founded background: the ordinary Cartesian 
partial derivatives d x , d y are replaced by the complex ones d z = (d x — id y )/2 and = (d x + id y ) / 2. 
In [[24l |25) three-dimensional derivative operators are introduced that behave similar to complex 
derivatives, that is, they are compliant with the rotation behavior of spherical harmonics in 3D. In 
[[26l[27l the Fourier transform of SE(3) is used in the context of engineering applications. For 
the efficient computation of SE (3) -convolution functions are expressed in terms of the unitary 
irreducible representations (UiR) of SE{3). The present work proposes a kind of intermediate rep- 
resentation. While the full SE (3) -UiR representation decomposes also the spatial variable in terms 
of Bessel-functions, the representation in terms of SO (3) -UiR leaves the spatial part untouched. 
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1.2 Preliminaries and Organization 

Most of the mathematical notations and conventions are adopted from [3 J regarding the geometry 
and parametrization of SE(3) and 50(3). Regarding the irreducible representations the notations 
are similar to E51 . 

In Section [2] we give a short introduction into the geometry of SE(3) and fix the related con- 
ventions. Section [3] introduces the necessary background of representation theory of SO (3). The 
main contribution of this work is presented in Section [4} where the left-invariant vector-fields of 
SE(3) are expressed in terms of the irreducible representations of SO (3). Finally, in Section [5] 
applications in the context of DWI and object detection are proposed. 



2 Life in SE(3) 

The motion group is the semidirect product of the rotation group 50(3) and the translation group. 
An element g e SE(3) is composed of a translation vector (x, j/,z) = rGK 3 and a rotation matrix 
R G SO (3) with multiplication law 

(r, R)(r', R) = (Rr' + r, RR) 

The Lie group SE(3) is generated by the six dimensional Lie algebra se{3) = TSE(3) spanned 
by the six left- invariant vector fields {T X) T y) T Z) J X) J y) J z } corresponding to the three transla- 
tions and the three rotation axis. They generate the right regular motion of smooth functions 
(j) : SE{3) 1 — y C, i.e. if we define (1Z h (fi)(g) := </)(gh) 9 then the application of a left-invariant 
vector-field A gives 

(Jl m 4>)(g) = ^ 

01 t=0 

where h : R H> SE(3) is a smooth curve with h(0) = e. And hence 

<f>{gh) = ll h (f)(g) = exp(A(g)t)<f>(g) 

Note that the tangents A{g) £ T g SE{3) do vary over the group and depend on the point g. The 
vector-fields are usually expressed in Euler angles. We parametrized the rotation in Euler-angles 
a, /?, 7 in ZYZ convention as follows 

Rg R^^R^^R^jQ 

where all rotations are counter-clockwise. In this parametrization the left-invariant vector-fields of 
the translation take the form: 

' % 

T y I = R*V 



(MM = i 



4(gh(t)) 



with V = (d X) d y) d z ) and for the rotation 



J x = cos a cot (3d a + sin ads — CQS a g 

sinp 

• 00 o sin a 

j v = — sm a cot po a + cos olOr H -cL 

sinp 
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Note that the fields depend explicitly on the position in the group. The conventions are the same to 
ones used in 0. 



3 Unitary Irreducible Representations of SO(3) 



It is well known that the spherical harmonics (see Appendix |7.1| for definition) are a orthogonal 
basis for the square integrable functions on the unit sphere L 2 (5 / 2 ). They have a variety of gentle 
properties with respect to rotations. The Wigner D-matrices play the same role for square integrable 
functions h 2 (SO(3)) on the rotation group SO(3) itself. They are representations of SO(3) and 
are written in Euler angles as 

(pi{g)) nm = 2^,(7, 0, a) = d^e-™ (1) 



where the real 'small' d-matrix d J nm is related to the Jacobi polynomials (see Appendix 7.3 ). As the 
D J (#) are representations of SO (3) they obey the multiplication law 

Ty(gh) = Ty(g)Ty(h) 

for any g,h e SO(3). For each order j G N they work on complex 2j + 1-dimensional vector 
spaces, i.e. — j < n, m < j. Furthermore, they are the irreducible representations of 50(3), i.e. 
there is no linear transformation A such that A T D J (g)A is block-diagonal for all g e S0(3). For 
j = 1 there is the direct relation to the original rotation matrix by 

/-I i 

D 1 ^) = SR^S T with S = — 0^2 

^ 2 V -1 -i 

Note that S is unitary, i.e. S T S = I. The irreducibility has an important consequence. By the Peter- 
Weyl theorem the irreducible representations form a complete orthogonal basis set of functions with 
respect to the group dot-product: 



/ f i 87T 2 

JqeSO(3) Z J t- 1 



r-2 

lgeSO(3) 

Thus, we can write any (j) e h 2 (SO(3)) as 



1 



8tt 2 

i=0 n=-j ro=-j 



where the expansion coefficients can be obtained by a simple projection 

L = (0,^nm>SO(3) = / <K#) D J nm (g) 

Jg£SO(3) 

onto the Wigner D-matrix. 
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3.1 Clebsch Gordan Coefficients of SO(3) 



The Clebsch Gordan (CG) coefficients interrelate the irreducible representations of different order. 
We denote a Clebsch Gordan coefficient by (Zto|ZiTOi, l 2 m 2 ) where /, /i, l 2 are different orders such 
that the triangle inequality \l\ — l 2 \ < I < \h + k\ holds, otherwise the coefficients vanish. The 
basic equation connecting two different representations of order li and l 2 is the following 

D Ln = Yl ^n^TO^Kl^l.^)^!/^!,^) (2) 

m^-\-m2 =m 
n l~l" n 2 =n 

Note the additional selection rule of the CG-coefficients, they only contribute if m = mi + m 2 . One 
also knows that (j0| jiO, j 2 0) = if j + ji + j 2 is odd. There is a variety of orthogonality and sym- 
metry relations for the Clebsch-Gordan coefficients making them itself a similarity transformation. 
We mention here the most important orthogonality relations: 

oo j 
3=0 m=-j 

^2 (jrn\j 1 m 1 J 2 m 2 )(fm f \j 1 m 1 J 2 m 2 ) = 5j^5 m ^ (4) 

777=7771+7772 

For example, suppose a sequence of numbers a J m with j < J is given. Then, for two fixed ji , j 2 
obeying \j x - j 2 \ < j < \j ± + j 2 \ we can compute f£{% 2 = Ej,m a{ n (jm\j 1 m 1 j2m 2 ), which 
contains all the information about the original a J m in an unitary way. Besides the orthogonality 
relations there are numerous other symmetry and associativity relations for the Clebsch Gordan 



coefficients, some of them are listed in Appendix 7.2 



3.2 Solid Harmonics and Spherical Derivatives 

In terms of the associated Legendre polynomials the components (Y J ) m = of the Racah- 
normalized spherical harmonics (for further details see Appendix |7.1[ ) are written as 

^(/? ) 7) = ^|^|^(cos(/5)) e i - 

Instead of /5, 7 defining a point on the sphere, we write in the following n £ S 2 as normalized 
Cartesian vector. There is a close relation between the Wigner D-matrix and the spherical harmon- 
ics: 

D^)Y^'(Rjn)=Y>) 

for any g G SO (3) and n <E S 2 . That is, the expansion coefficients of a spherical harmonic 
expansion rotate by the application of Wigner D-matrices. With this we can identify the central 
column of a Wigner D-matrix with a conjugate spherical harmonic by rotating a spherical harmonic 
along the z-axis 

iy(g)Yt(e z )=Yt(R g e z ) 
and with the additional knowledge that (Y j (e z )) m = 8 m ^ we have 

Dy 7 ,/3,0)=^(/3, 7 )- 



5 



Next to the spherical harmonics, the so called solid harmonics 



Rl(r) = |if Y^(r/|r|) 

are solutions of the homogeneous Laplace equation, i.e. AR J m = 0. They are homogeneous 
polynomials of degree j, that is, R J m (Xr) = X j R J m (r) for any AgR. So, we can define a differential 
operator as follows 

dl := /&(V) 

which is a spherical tensor operator, i.e. it inherits all rotation properties of irreducible representa- 
tions. Further note that d 1 = S V. 

4 Diffusion Equations in terms of the Irreducible Representa- 
tions of SO(3) 

The quadratic forms in the left-invariant vector fields [3 J generate the 'dynamics' on functions 
cj) : SE(3) \-> C. Our goal is to find them in terms of the irreducible representations of 50(3). 
More precisely, suppose we have an evolution equation 

dt<Kg,t) = H</>(g,t) 

where g is parametrized as proposed above and the evolution generator is a polynomial in the left- 
invariant vector fields H = H{T, J) = H(T X) T y) T Z) J X) J y) J z ). That is, H acts as a differential 
operator in d X) d y) d z and <9 a , c^, d 1 on the function (p. By decomposing cp in terms of Wigner 
D-matrices 

1 J 

3=0 n=-j m=-j 

we will be able to show that the evolution equation in terms of /^ m (r, t) is just a differential operator 
in the spatial coordinates and algebraic for the angular coordinates. The form of the equation will 
not depend on the particular choice of the chart chosen to parametrize SO (3). That is, we expect 
the equation to be 

where the H 3 ^ ml are differential operators in the spatial coordinates. To find H we have to express 
all appearing quantities in terms of the irreducible representations D J nm , i.e. we have to evaluate the 

matrix elements (HD J nm , D J n , m ,) S o(3) '•= H 3 ™ m , or (Hep, D J nm ) S o(3), respectively, which is subject 
of the next section. Note, that by the help of the UiR of the full group SE(3) [26], the above equa- 
tion can be made purly algebraic as long as no external gauge field is used. However, this approach 
would also approximate the translation/spatial part of the function, which is in our context not nec- 
cessary and avoids any approximation artefacts coming from this side. Of course, the Wigner-D 
expansion is restricted to a finite cutoff index j < J in practice, and thus, we also have approxima- 
tion errors. More precisely, if Pj is the linear orthogonal projector onto the subspace of functions 
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spanned by the Wigner-D matrices up to an order of j < J, then the approximated evolution gen- 
erator is Hj = PjHPj and the corresponding evolution operator is exp(Hjt) = exp(PjHPjt). 
That is, the obtained time evolution is not the projection Pj exp(Ht)Pj of the original one, which 
is important to note. 



4.1 The Left-invariant Vector Fields 

We start with switching to the complex representation which is obstructed by the irreducible repre- 
sentation of rank j = 1: 

/T_i\ f -(%-!%)/ V2\ 

T= To = % = SRjV = D 1 (y) T 1 (5) 

\r + ij \-(% + iTy)/V2 J 




-(Jx-'lJy)/V2 

J = | Jo I = I I , (6) 

-(X + iJy)/V^ 

where we formally discriminate between both representations by numeric or character x, y, z subindices. 
The J±i are well known from quantum mechanics and engineering literature ll27ll28ll29l[30l[3Tll32ll 
as the body-fixed rigid rotor angular momentum operators. They obey the following equations: 

&DL = -™z>L (7) 

J±iDL = VjO" + l)/2-m(m±l)/2Di (m±1) . (8) 

Thus, we can directly compute the action of the operators J^onto the coefficients fields /^ m (r). 
Therefore, we denote the collection of all coefficients by a bold Latin letter f and we access ele- 
ments of f by round brackets, i.e. (f) 3 nm = fl m and thus, the coefficient {J z i ) J nm is the projection 
of {J z 4>) onto the orthogonal Wigner D-matrix basis: 

WYnm = (Jz(p,D J nm )so(3) = / dg R (Jz<I>){9t9b) D J nm {g R ). 

JSO(3) 

where the integration ranges over SO (3). With formula ([7]) and partial integration we can find 

(JMm = ~ j dg R m J»BU9r) = im ft n (9) 
and similarly with equation ([8]) we can proceed as follows: 

( J±if y nm = ~ J dg R 4>{g) J±iD J nm (g R ) 

= -Vi(i + l)/2-m(m±l)/2 J dg R <f>(g) D j n{m±l) (g R ) 

= + l)/2 - m(m ± l)/2 f n{m±l) (10) 
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The translations are a bit more intricate to compute. We apply the right hand side of equation ([5]), 
i.e. 7 = T> l {g) T d 1 onto the field (j) and project onto the orthogonal Wigner basis: 

W) 3 nm = J ' dg R (T k (f))(g) D J nm (g R ) 

q=-l,0,l J 



By using the integral formula (49) for triple products of Wigner D-matrices we get 



(WL = 8^ E ^' + 1 ) d lfL'Jd9RDl k Di m/ Di m 

j' ,n' ,ra / ,q 

= E E (ID 

j f =j — l,j,j + l n—n' + q 
q— — 1,0,1 m=m / + fe 

which is already our first main result. It gives the action of the translational left-invariant vector- 
fields Tk in terms of the Wigner expansion coefficients fl m . Note that the sums are very sparse, 
which is from a computational viewpoint quite important. The sum for / runs only over three 
terms. 

4.2 Quadratic Forms 

While the linear forms in the left-invariant vector-fields generate Euclidean motion, the quadratic 
forms generate diffusion. The important Laplace Beltrami Operator on SO(3) 

J 2 = Jl + Jy + Jl = J+lJ-l + JO + J-l J+l (12) 

generates diffusion on SO (3). The action of J 2 onto the coefficients fl m can be computed with 
the help of equation ([9]) and ([10]) to 

(J 2 f) J nm = -jti + l) fnm- (13) 



Other arbitrary products Jk'Jk can be computed quite easily by the use of equation ([9]) and ( [TO] ). 



Also products of rotations J k and translations %* can be computed by the use of the first-order 
equations ([7]), ([8]), ( pT) . You can find the detailed formulas in Appendix |7.4[ However, the product 
of two translations is again more cumbersome to evaluate. The vector field %t Tk transforms with 



respect to k and k' like a Cartesian rank 2 tensor. By using equation ( |32| ) we can write it in terms 
of spherical tensors as follows (for a proof see Appendix |7.5[ ): 



%% = ^ - ^ E ( lk '\ 2 ^ lk ) (( d2 ) T ^ 2 )p ( 14 ) 



3 3 

v- 



where A = d 2 x + dy + d 2 z , which is related to the trace of Tk'Tu, i.e. A = Ylk \T~k\ 2 - The second 
term is related to traceless matrix TyTk — A/35 kk *. By using equation ( [T4| ) we can proceed like in 
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the linear case: 

(%%fy nm = I dg R (7k'Tk<f>)(g) D 3 nm (g R ) 



g,p=-2,...,2 



= f/L"^ E ^^<l^|2p,lfe>0n|jV J 2g>0m|jW,2p>^/i; m , 



k f =p-\-k 
j>=j-2,...,j + 2 
n—n 1 -\-q,m=m / -\-p 



where we have again used formula ( [49] ) for the triple products of Wigner D-matrices. To simplify 
the formula we consider the special case k = k f — 0, i.e. |7o| 2 = 7^, which gives 

(7?f)L = + | E |y^<^|jV,2g)(jm|j'm,20)^/^ m . (15) 

j / =j-2,...,j+2 
n—n f -\-q 

On the other hand |71i| 2 + |7I| 2 = 7^ + 7^ gives 

(C£ + 7?)f)L = ^/L-f E ^^(^bV,2g}(jm|/m,20}a^; m .(16) 

j'=j-2,...,j+2 

Note that there is only a small change between T£ + Ty and 7^. One can easily see that + 7^ 2 + 

r 2 = A. 



4.3 Restriction to R 3 x S 2 

Given a function (j) : IR 3 x S 2 ^ C it is natural to extend to functions <\> : SE(3) \-> C by 

0(r,R) = 0(r,Re,) = 0(r,n) 

Hence, there is the symmetry 0(r, R z?a ) = 0(r, I) for any rotation around the z-axis. In Euler 
angles this means that 0(r, 7, /?, a) = 0(r, 7, /3, 0). By this restriction the evolution operator H has 
to be invariant with respect to a rotation around the z-axis 

H(T, J) = H(R Zja T, H Z ^J) 

(see for details). Note that such an invariant H always depends on T Z) T£ + and J% + Jy 
(at least if no other external gauge image is used). To understand the implications for the Wigner 
D-expansion we have to look for Wigner D-matrices that satisfy £^ m (7, /?, 0) = D J nm (j, /3, a) for 
any a. From the explicit form in equation ([T]) it is easy to see that this holds for all D 3 nQ . So, the 
expansion in terms of Wigner D-matrices reduces to the usual spherical harmonic expansion 

1 J 

<Kr, 7 ,/3) = 0(r,7,/5,O) = ^EE E ^' + *) & °) A 

3=0 n=-j m=-j 

1 i 

J=0 n=-j 
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To translate the general action of the translational vector-fields in equation ( pT) onto this special 
case, we just have to set m = m! = 0. The selection rules of the Clebsch-Gordan coefficients imply 
k = 0, meaning that only the % gives non-vanishing contributions, which is in agreement with 0) 
that only those operators survive the construction that are well defined on the cosets SO(3)/SO(2). 
For this special case we visualize the matrix elements graphically in Figure |4~3j where we use the 
abbreviation 



Z J Sj ,, m , = (Z 3 jf ) nn , = Kn> +q %^{3n\j'n\ Jq)(j0\j% JO) d J q 

q =-j,...,j ZJ + 1 



for the block matrices on the upper and lower secondary block diagonal, where J is the differential 
order of the operator. The rank j is the rank of spherical tensor field which is returned after appli- 
cation of 7a 3 --, and f is the rank of the incoming spherical tensor field. So, we can write equation 
( |TT| ) for k = more compactly as 

J+l 3-1 

(7()f)n = ^){j+l),nn' + ^ ^]{j-l),nn' f n > 

n'=-{j+l) n >=-{j-l) 

= (zW i+1 + z iW i_1 )„ ( 17 ) 

Note that the matrices Z are also very sparse, the non-zero elements are shaded in blue in Figure 



43} In the same fashion we can write the second order operators 

A • 2 
-P + - 
3 3 



(W = ^F + ^Z^f^ + Z^.F + Z^F- 2 ) (18) 



and 

(YT 2 + T^fV = - . 

3 3 



(Ctf + W = ^ fi -?( Z ,V 2 ) fi+2 + Z ?/ i + Z ?0--2)f i - 2 ) d9) 
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Both operators generate horizontal and orthogonal diffusion ([3 J). There are also generalizations 
of these diffusion generators by convolutions with rotation symmetric functions on £2. These 
convolutions are known to have a very simple diagonal form in terms of spherical harmonics. If we 
denote the convolution operator by C, then the application is defined to be (Cf ) J = cjf J , where Cj 
are the expansion coefficients of the rotation symmetric convolution kernel on 5 2 . Two operators, 
which can be derived from a Tikhononv-regularization problem are C T 7^C and on the other hand 
T Z C T CT Z . They can easily expressed in our framework: 

(C T T 2 C f)' = ai jP + 2 -{a 2 Z 2 [j+2) P +2 + aaZ^-P + a,Z 2 {j _ 2 f~ 2 ) (20) 

where the factors are a\ = a 3 = \cj\ 2 ,a 2 = c]c j+2 and On the other we have 

9A 2 

(T z C T C%i) j = a,—P + -{a 2 Z% +2) P +2 + a % Z%P + a 4 Z 2 j(j _ 2 f- 2 ) (21) 

where a\ = a 3 = (|cj + i| 2 + |cj_i| 2 ), a 2 = |cj + i| 2 and a 4 = |cj-i| 2 . Note that the same operators 
can also be derived for orthogonal diffusion (T x 2 + 7^). 

5 Applications and Experiments 

We want to show two examples where the introduced framework can be applied. First, an exam- 
ple in MR-imaging will show that the introduced operators can effectively used as regularization 
terms for the estimation of fiber orientation distributions on the basis of HARDI-imaging. In this 
experiment the SH-based approach will be compared with its discrete counterpart. And secondly, 
we use the translation operator to implement a 3D circular Hough transform efficiently. But before 
starting, some details on the discrete implementation of the derivative operators are given and some 
general issues are discussed. 

5.1 Discrete Implementation 

All differential operators are implemented by ordinary finite difference schemes. To keep the 
implementation efficient we restricted our implementation to first order approximations. For the 
discretization of the first-order derivative operators Zj ., we use a simple central finite difference 
approximation, that is, the convolution kernel of a partial derivative along an arbitrary axis reads 
[— 1 l]/2. We also explicitly implemented second order derivatives T?--, by the common difference 
scheme: For example, the convolution kernels of d xx and d xy are 














) — 


1 





-1 


d X x 


1 


-2 


1 






















-1 





1 



the others can be obtained by permutation and symmetry. From these rather rough approximations 
we cannot expect too much. In particular, the forward Euler-integration of first-order central dif- 
ferences is usually a no-go. Unfortunately, usual solution for convection dominated problems by 
upwind/downwind schemes cannot be applied due to the global nature of the harmonic represen- 
tations. Nevertheless, we want to rely on these crude approximations for the sake practicability. 



11 




position (voxel) position (voxel) 



Figure 1: The kernel 7o was iterated 150 times with a step- width of 0.05 resulting in a translation 
of 7.5 voxels. The initial condition was set to 0o(r, n) = ^ 2 S Uz (n). On the left the maximum 
max n (p(zn zi n) along the z-axis is plotted. On the right the mean-component j = 0, i.e. f°(zn z ) 9 
is visualized. 

Of course, there are several more sophisticated approximation schemes (J33l [34l [35]]), but due to 
the inherent slow-down of the computation speed, in particular in 3D, we only consider the most 
simple scheme. And we will see in the experiments that it is enough to obtain reasonable good 
results. 

In a first small experiment, we investigate the behavior of the simple forward integration of % 
in detail. Additionally, we look at the forward integration of 7o + 0.17^ 2 , i.e. a small diffusion 
term is added to make the results more stable. In particular, we used o (r,n) = / 2 8 Uz (n) 
as initial condition and iterated the integration 150 times with a step- width of 0.05 resulting in a 
translation of t = 150 • 0.05 = 7.5 voxels. To analyze the result we have a look at the maximum 
of the orientation distribution along the z-axis, namely max n (j)(zn z , n) and the j = component 
f°(zn z ) 9 which is the mean along the orientation coordinate. For reference, we made a simple ID 
experiment without any orientation involved. It was just integrated (j)(z) \-t + 0.05 • d z (f)(z) 
with the same first-order central difference approximation of d z . In Figure[T]the results for different 
cutoffs L > j and the ID experiment are shown. The ID reference experiment shows the typical 
oscillations. Similar oscillations can also be observed for j = component. In particular, for low 
cutoffs the oscillations are pretty heavy and converge for large L towards the ID reference. For the 
maximum of the orientation distribution the behavior is different. Instead of oscillations a tail is left 
over, while certain peaks appear where the oscillations of j = component dominate. Looking at 
Figure [2] where the small additional diffusion term is considered, one can observe less oscillations 
and a more regular behavior. In Figure [3] we show the resulting orientation fields for L = 12 in 
glyph representation with and without diffusion. Once again one can see the spike for the version 
without diffusion and the more regular behavior for the diffusion-regularized version. 

5.2 Spatially Regularized Spherical Deconvolution 

Magnetic resonance imaging (MRI) has the potential to visualize non-invasively the fibrous struc- 
ture of the human brain white matter ll36l . Based on the directional dependency of water diffu- 
sivity in fibrous tissue it is possible to reveal underlying connectivity information. The accurate 
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Figure 2: The kernel 7o was iterated 150 times with a step-width of 0.05 resulting in a translation 
of 7.5 voxels. The initial condition was set to 0o(r, n) = e"l r l ^ 2 S Uz (n). On the left the maximum 
max n (j)(zn zi n) along the z-axis is plotted. On the right the mean-component j = 0, i.e. f°(zn z ) 9 
is visualized. 





Figure 3: The same setting like inland [2] but the orientation distributions are visualized by glyphs. 
The spherical harmonic cutoff is set to L = 12. The underlying gray- values indicate the maximum 
of the orientation distribution. 
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and reliable processing and estimation of fiber orientation distributions is a major prerequisite for 
the processing of such data. There are numerous methods for estimating orientation distributions 
on the basis of the diffusion-weighted MR-signal. We will focus spherical deconvolution Il37lh 
which is one way to estimate the so called fiber orientation distributions (FOD) on the basis of the 
diffusion-weighted MR-signal. The idea is based on a model-driven deconvolution scheme to turn 
the diffusion weighted MR-signal into a FOD. The goal is to find a FOD / such that 



J(f) = J J |(H/)(x,n) - S(x,n)\< d*dn (23) 

R 3 xS 2 

is minimized. Here S denotes a quantity derived from the MR-measurement, e.g. in standard q- 
ball imaging just the ratio M(x, n) /M (x), where M is the diffusion weighted image and M the 
measurement without diffusion weighting. The operator H denotes the spherical convolution with 
the so called fiber response function: 

(H/)(x,n)= / /i(n.n , )/(x,n , )dn , J (24) 

Js 2 

The fiber response function is typically chosen to be a function of the form h{t) = exp(— At 2 ), or 
is estimated from the measurement itself. The main problem is that H is practically not invertible. 
One way to solve the problem is to introduce a non-negativity constraint [10], but in case of low- 
quality data the problem is still hard to invert. Therefore, regularization techniques JH |3l [El were 
introduced, which make the solution unique and stable. An additional term R(f) is added to the 
cost function as: 

J Kg (f) = J(f) + R(f) 

The fiber continuity/contour enhancement kernel provides such an additional cost function that 
perfectly fits to the nature of fiber orientation distributions: 

R FC (f) = ^JJ(^' V/) 2 rfxrfn. = A J J (To/) 2 rfxrfn. (25) 

R 3 xS 2 R 3 xS 2 

We can expect from this kernel that it prevents 'arbitrary' smoothing and that it will preserve and 
emphasize the fibrous nature of the data. As the regularizer includes the spatial neighborhood of 
each voxel, we have to be careful at the transition area between gray and white matter: a boundary 
condition is needed. Instead of a hard boundary condition we decided to keep the FOD small in 
the background area (non white matter). Thus, we have an additional term in the cost function that 
suppresses the FOD in the background, 



^mask(/) = A mask J J /(x, n) (1 - w m (x)) dxdn, (26) 

R 3 xS 2 

where w m (x) is a white matter mask. That is A mask acts as a suppression strength A mask for gray- 
matter. We found a value of 1 to be a good value. 

The choice of regularization strength is a crucial issue. We found that too strong regularization 
emphasizes discretization artifacts, which is shown by a violation of the rotation invariance, that 
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is, the results depend on the absolute directions of the bundles, with an unpredictable dependence 
on the sphere discretization and the underlying voxel grid. On the other hand too low values lead 
to less stable results. We found by a simple visual inspection of a simulated crossing a value of 
A = 0.005 to be a good trade-off. 

5.2.1 Optimization and Implementation 

In order to find the optimum of J reg , we have to compute the variation of the objective and set it to 
zero for the necessary condition for the optimum: 

5 j J reg = H T H/ - H T S - \T*f - A mask (l - w m )f = (27) 

which leads to 

(H T H - AT 2 - A mask (l - w m ))f = H T S (28) 

To solve this linear equation we employed an ordinary conjugate gradients scheme. The implemen- 
tation of the operator on the left-hand side in terms of spherical harmonics is based on equation ( [T8] ). 
The operator H is also quite simple to implement in terms of spherical harmonics, because it is a 
diagonal matrix. The elements on the diagonal are just h(t)Pj(t)dt. If we want to implement 
this with a discretized sphere the operator H has to approximated by interpolation. We just used 
spherical harmonics to do this interpolation, in the same way like in [35) the spherical Laplace- 
Beltrami operator was approximated. The operator Tq is discretized very similar to our spherical 
harmonics implementation with finite differences (comparable to ll35lh Appendix F). Throughout 
this experiment the sphere was discretized with 512 direction, which were determined such that 
the energy of the configuration of 512 electrons on the surface of a sphere is minimized, where the 
electrons repel each other with a force given by Coulomb's law. 

For both methods the conjugate gradients algorithm was iterated 100 times, which was enough 
for convergence. 

5.2.2 Experimental Setup 

The goal of the experiment is to compare the proposed spherical harmonic representations with the 
discrete angular representation of the fiber orientation distributions. Therefore, we simulated the 
MR-signal with 64 gradient directions of a crossing region. The signal was simulated by using the 
standard exponential model £ nfib (n) = e - bD ( n - n ^) ? that is, no diffusion perpendicular to a fiber 
is assumed. We have chosen bD = 1, which emulates a b-value of 1000s /mm 2 and a typical 
diffusion coefficient for the human brain. The generated signal was distorted by Rician noise 

S'noisy = ^(S + n rea i) 2 + ^! 2 ma g, where n rea i and n imag are normally distributed real numbers with 

standard deviation a. The signal-to-noise ratio is defined as SNR = 1/a, that is, the SNR is 
calculated with respect to the b = measurement. The crossing was created on a 24 x 24 voxel 
grid, where the tracts of the crossing are on average 5 voxels thick. To get an impression look at 
Figure [4j 

To measure the performance of the deconvolution method, the local maxima of the estimated 
FODs are extracted and compared to the ground truth direction. To efficiently search for a lo- 
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Figure 4: Two simulated crossing situations to compare SH representation with the angular discrete 
version. The simulation was performed with 64 gradient directions at a b-value of 1000s /mm 2 and 
with a diffusion coefficient of 10~ 3 mm 2 /s 



16 



cal maximum we followed the this procedure: We made a Voronoi tessellation of the sphere[]and 
search for those direction whose value is above all its neighbor with respect to the tessellation. To 
accurately determine the underlying continuous direction, a quadratic form is fitted to the neigh- 
borhood of the putative local maxima and the maximum of this quadratic form is used as detected 
direction. A ground truth direction is said to detected, when it is in a range of 10° from a detected 
direction. To measure the performance we used precision, recall and the f-score0 To generate the 
performance measures the simulations were repeated 100 times. 

5.2.3 Method Comparison and Discussion 

In Figure [4] we compare visually the spherical harmonic implementation for a cutoff of L = 8 
with the discrete version with 512 directions on the sphere. We consider a crossing angle of 35° 
and 45°. For both situations one fiber bundle direction was chosen along the underlying Cartesian 
coordinate axis. For the larger crossing, which is easier to resolve, we assumed a relatively low 
SNR of 7. For this case one can see that the discrete version is more susceptible to noise than the 
spherical harmonic version, which is not astonishing due to the implicit regularization by the finite 
SH-cutoff of L = 8. For the smaller crossing angle of 35° a higher SNR of 20 was assumed. In this 
case a SH-representation of L = 8 is nearly at its limits to discriminate between both directions, 
while the discrete version can still well distinguish. Further, one can observe that for both methods 
the FODs along the horizontal Cartesian axis are sharper than along the skew axis. But this effect 
is more prominent for the discretized version than for the spherical harmonic representation. 

In Figure [5] we show quantitative results. The crossing was simulated for varying crossing 
angles between 30° and 90°. Additionally we varied the absolute pose a of the crossing. For 
a = 0° the horizontal bundle is aligned with underlying x-axis of the Cartesian grid. With growing 
a the whole configuration is rotated clockwise. The crossing was simulated at a SNR of 50, which 
is a realistic scenario. As a baseline experiment we show results of the so called Constrained 
Spherical Deconvolution (CSD) approach [Toll , where an additional positivity constraint is used 
to obtain more stable results. Apart from the unregularized CSD approach the results obviously 
depend on the absolute angle of the configuration. The discrete approach is able to cope quite good 
with small crossing angles, but shows independent of the crossing angle less precision than the 
SH-based approach. In particular, for crossing angles above 45° the SH-based approach solves the 
task nearly without any error. 

In Figure [6] the results are further investigated by scatter plots in the 0, #-plane. With the same 
parameters as above a crossing of 50° was simulated and reconstructed by the three methods. Each 
detected direction is indicated by a small dot at its corresponding angle 0, 6. The crossing was 
simulated twice, for an absolute angle of a = 0, i.e. one direction is along the x-axis, and secondly 
for an absolute angle of a = 15°. For a = the SH-based approach is able to resolve the direction 
along the x-axis (cj) = 90°) perfectly, while the other direction (cj) = 40°) is a bit more blurry. On the 
other hand, the discrete approach has severe problems with the cj) = 90° direction, which explains 
the lack of precision, i.e. apart from the true direction there are some additional local maxima that 

Un the case of the SH-representation the orientation distribution was sampled with the same directions like the 
angular discrete algorithm uses 

2 Let TP be the number of successfully found ground truth directions, let FP be the number of detections that are not 
in a range of 10 degree to a ground truth direction, and let FN the number of ground truth direction that are not detected, 
then precision = TP/ (TP + FP) and recall = TP/ (TP + FN) and f-score = 2 precision • recall/ (precision + recall). 
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Figure 5: Detection accuracies in terms of precision/recall and f-score for different crossing angles 
and absolute angles a. The measurement was simulated at a SNR = 50 with a value of bD = 1 
and 64 gradient directions. 
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Figure 6: Scatter plots in the (f)-9 plane for the crossing configuration with an crossing angle of 50°. 
The angle cj) is plotted along the x-axis, 6 along the y-axis of the scatter plots. The intersections of 
the thick red lines indicate the expected ground truth directions. The dotted lines indicate the 10° 
detection tolerance. On the left, the results for the configuration with an absolute angle of a = 
are shown, on the right the results for a = 15° are given. 

produce false positives. The main reason is the interplay of the 64 gradient directions and the 512 
discrete directions of the FOD. The effect is reduced by an increase of measurement directions. 
But also the SH-based approach has problems when the number of measurement directions is too 
low. They are revealed for an absolute angle of a = 15°. Besides the uncertainty caused by 
the measurement noise one can observe a systematic bias. For example, for the direction along 
(f) = 55°, the center of the distribution is shifted in 9 by approximately 5°. Also for the other 
direction the distribution is a bit squeezed. We also found that the main reason is the low number 
of measurement directions. For example, for 128 gradient directions the estimated directions are 
unbiased. Another way to reduce the effect is to decrease the expansion cutoff of the spherical 
harmonic representation. For L = 6 and 64 gradient directions the estimates do not show a bias. To 
conclude the differences: both methods have problems when the number of measurements become 
too low. While the discrete approach shows scattered, multimodal distributions, the SH-based 
approach shows a slight systematic bias but the distributions stay unimodal. 

In Figure [7] we show a real world example of the human brain. The setting of the measurement 
is nearly the same like in the simulations. A b-value of 1000s /mm 2 and 61 gradient directions were 
used with an isotropic resolution of 2mm. A kernel of the form S nnh (n) = e ~ bD ( n - n ^) w \fa frjj _ \ 
was used as a model for deconvolution. Overall, both methods work very similar. Figure [7] shows 
a coronal section in glyph representation. One can observe that the SH-based approach produces 
a bit more negative values, which are indicated in black. The green rectangle highlights a region 
where the differences are largest. The red rectangle shows a regions where the discrete approach 
shows a direction which does not appear for the SH-based approach. Whether the direction is true 
or not is difficult to say, but the fact that it precisely points along the x-axis makes it dubious. 



19 





rxi r 



//if 1 

* H \ \ t i f — 




77 

S \ \ H « * j 




Figure 7: A real world example of the human brain. A coronal section is shown. On top the results 
of the SH-based approach with L = 8, on the bottom the discrete approach with 512 directions. 
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5.2.4 Memory Consumption and Running Time 

The memory consumption of the SH-based and discrete approach is easy to compare. We want 
to consider the above real world experiment as an example. The whole volume has a size of 
96 x 96 x 60. As the FODs are symmetric we need to store only the even expansion coefficients 
of the spherical harmonic representation. Additionally, the FODs are real, hence we have to store 
j + 1 numbers per expansion field, instead of 2j + 1 like for a general complex field due to the 

symmetry Ym = (— l) m Yl m . Thus, each voxels consumes (L + 2) 2 /4 complex numbers, resulting 
in 50 • 8 bytes for L = 8 and double precision. Overall, one volume needs 96 2 • 60 • 50 • 8 = 220 
MB in SH-representation. On the other hand, in discrete representation with 512 direction needs 
96 2 • 60 • 512 • 8 = 2264 MB, which is 10-times more compared to the SH-representation. Recall, 
that the conjugate gradient algorithm needs four instances of the volume in memory. 

The running time is more difficult to compare. On the one hand, we have to compute the ap- 
plication of % which basically consists of finite differences, on the other hand, the convolution 
operator H has to be implemented. Let us consider the computation of % first. In case of the 
discrete approach one has to compute for each of the 512 components six second order finite differ- 
ences, which are linearly combined with weights depending on the corresponding direction. For the 
SH-based approach one also have to compute six finite differences, but the linear combinations of 
them to obtain the final values are more more expensive. In particular, we implemented separated 
functions for the operators Z 2 J+2 , Z 2 J _ 2 and Z| thus several values are computed repeatedly. 
In practice we found that one application of % with 256 discrete directions is comparable to an 
SH-based application of % with L = 10. 

The computation of the operator H is negligible in SH-representation, while it is the bottleneck 
for the discrete approach. Here, the running time heavily depends on the algorithm used for the 
matrix multiplication. In fact, the execution times can differ about a factor of 10. While the highly 
optimized BLAS matrix-multiplication shipped with MATLAB is quite competitive, a standard 
non-optimized version can slow down the running time dramatically. To give an example, to re- 
construct the above volume with 100 CG iteration with 512 directions takes on a Intel Xeon X7560 
@ 2.27GHz about 20 minutes with a highly optimized multiplication. On the other hand, with a 
standard BLAS implementation it needs above an hour. For comparison, our implementation for 
L = 8 takes about 10 minutes on the same machine. 

5.3 A Spherical Hough Transform 

The detection of spherical objects is an everyday issue in biological image processing. The Hough 
transform ll38l is often the method of choice. In this section we want to show how the contour 
completion kernel can be used to implement a variant of this approach. Suppose, we have a 
volumetric image of a solid spherical object, such that the gradient field g(r) at the surface of the 
object points away from the object center r . Let m(r) = |g(r)| be the gradient magnitude and 
v = g/m the gradient direction. If the object has radius p, then we know that for all n e S 2 
approximately v(r + pn) = n holds. This fact can be easily used to get a 'evidence' or 'voting' 
map for the center r of the object. Therefore, let 5 n (m) an indicator function on the two-sphere, 
not necessarily a delta- function but a bit blurred. Then, we count for each putative object center r 
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Figure 8: A toy example for the spherical Hough transform. A distorted shell was rendered into a 
volume. The top row shows the initial data. In the two bottom rows the initial gradient orientation 
distribution is integrated with the kernel 7o+0. 1 -7^ 2 , with a step width of Ap = 0. 1. The orientation 
distribution is expanded up to L = 4. The gray value background shows the j = component. 
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Figure 9: The evolution of the Hough transform for a real world example of confocal laser scan of 
airborne pollen grain. All images are maximum intensity projections of the 3D volume image. The 
Hough transform was applied with the same setting like in the toy example. 



how often v(r + pn) = n is approximately fulfilled for each possible direction nG5 2 : 

h(r , p) = m(r + pn) 5 n (v(r + pn)) dn, 

where each contribution is weighted by the gradient magnitude of its originating voxel. The local 
maxima of the map h give evidence for the presence of a spherical object with center r and radius 
p. Now, the integrand of the above equation can be expressed in terms of the horizontal translation 
operator: 

H(r, n, p) = m(r + pn) 5 n (v(r + pn)) 

= e~ pn - v (m(r) 5„(v(r))) = e~ pTo (m(r) 5„(v(r))) 

That is, we have a quite simple algorithm to get the Hough voting map: we initialize with if (r, n, 0) := 
m(r)5 n (v(r)) and successively Euler-integrate 

H(t, n, p( n+1 >) = if (r, n, p (n) ) + Ap AH(r, n, p(»>), 

where .4 = — %. To get a more stable response we added, as discussed above, a slight amount 
of diffusion A = — % + 0.1 • T^. Note that the proposed Hough transform only works for solid 
objects with surface gradients pointing away from the center. However, it is easy to switch to 
inward pointing gradients by using A = 7o + 0.1 • instead. 

In Figure [8] we show a toy example for a shell, that is, the voxel on the surface of the sphere 
are set to one. Hence, both gradients are present: inward pointing gradients at the outer border and 
outward pointing gradients at the inner border. We decided to let the outward pointing gradients 
to be translated towards the center, i.e. A = 7o + 0.1 • 7?. The shell in Figure [8] was rendered 
with a radius of 7 in a 32 3 -voxel cube. Every second voxel on the surface was deleted randomly 
and Gaussian random noise with standard deviation of 0.3 was added. To get an impression Figure 
[8] shows an isosurface at gray value level 0.6, a central slice and a maximum intensity projection 
(MIP) of the toy example. The orientation field H is expanded up to L = 4 and integrated with a 
Ap = 0. 1. The initial gradient field g was computed on a Gaussian smoothed image of width a = 1. 
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The lower two rows in Figure [8] show the evolution of orientation field H for p = 0, 2, 4, 6, 8, 10 
in glyph representation. The underlying gray value image is the final voting map h, i.e. the £ = 
component of the orientation field H. One can see how the gradients at the outer border are shifted 
outwards and gradients of the inner border are translated towards the center. Approximately for 
p = 7 the voting maps looks as desired. 

6 Conclusion 

This article worked out the formulation left-invariant convection/diffusion equations on SE(3) in 
terms of the irreducible representations of SO (3). From a computational viewpoint the main ad- 
vantage of the proposed formulation is the low memory consumption in comparison to a angular 
discretization of the two-sphere or rotation group, respectively. With a low number of basis func- 
tions the functions are well described without hurting the rotation covariance. The DWI example 
showed that even with 10 times less memory consumption the results are still comparable to the 
angular discrete approach. In terms of accuracy both approaches have their own problems. 

Applications to the full group SE{3) remain subject to future work. For example, the detection 
of helical structures in cryo electron micro-graphs ll39l might be a good playground. We further 
plan to provide the elementary operators as a open source toolbox to give the scientific community 
a chance to try the proposed framework with a small amount of effort. 

7 Appendix 

7.1 Spherical Harmonics 

We always use Racah-normalized spherical harmonics such that Y^(r) T Y^(r) = 1, or Y^(r) T Y^(r / ) = 
i^(cos(r, r')), where the Pg are the Legendre polynomials: 

W) = ^,af(t 2 - 1)'. 

In terms of the associated Legendre polynomials the components of the spherical harmonics are 
written as 

YttM) = y / f^ p r(cos(0))e i ^ 
Mostly we write r G S 2 instead of (0, 9). The Racah-normalized solid harmonics can be written as 

Rl(r) = ^ + m)!(^-m)!^ ^^ 2 ^ (x - i y y(-x - iy)^, 

where r = (x, y, z). They are related to spherical harmonics by R £ m (r)/r £ = Y^(r) 
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7.2 Clebsch Gordan Coeffcients 

The Clebsch Gprdan coefficients of SO(3) fulfill several orthogonality relations: 

^2{jm\j 1 m 1 , j 2 m 2 ) (jm^m 1 ^ j 2 m' 2 ) 

^mi ,777^ ^7772 ,777-2 

j,m 

2j + l 



^7771 ,777^ ^7772 ,777 2 
j,777 J1 

^2 (j m \hmij2rn2)(fm f \j 1 m 1 J 2 rri2) = 5 jd >5 m , m > 

777=7771+7772 

^U m \ji m iJ2rn 2 ){jm\jim 1 J' 2 'mi 2 ) = |i^-^ 2 j^m 2 , 

mi,m J2 



For particular combinations there are simple formulas 

(£m|(£- X)(m- fi),Xfi) = 



e + m\ 1/2 f i-m\ 1/2 f 2i^~ 1/2 



A + /i I y A — n J y 2A 
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There are several symmetry relations 



{jm\j 1 m 1 ,j 2 m 2 ) = (jim 1 ,j 2 m 2 \jm) 

{jm\jim u j 2 m 2 ) = {-l) 3+31+n (jm\j 2 m 2 , jimi) 

{jm\j 1 m 1 ,j 2 m 2 ) = {-l) 3+n+n {j{-m)\j 1 {-m 1 ), j 2 (-m 2 )) 

(jm|jimi, j 2 m 2 ) = J ^ * * ( - 1 ) J1 +mi ( j 2 m 2 1 jm , ji ( - m 1 ) ) , 
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and associativity relations: 



(J,M\j 1 +j 2 ,m 1 + m 2 ,j 3 ,m 3 ){j 1 + j 2 ,m 1 + m 2 \ji,m 1 , j 2 ,m 2 ) = 

{J,M\j 1 +j 3 ,m 1 +m 3 ,j 2 ,m 2 )(j 1 + j 3 , mi + ra 3 |ji, mi, j 3 , m 3 ) (43) 

where J = ji + j 2 + j 3 and M = mi + m 2 + m 3 . And for j 3 > ji + j 2 we have another one: 

(j3 - Ji — J2, mi +m 2 + m 3 |J - ji,mi + m 3 , j 2 ,m 2 ){j 3 - j 1 ,m 1 + m 3 \j 1 ,m 1 , J,m 3 ) = 
(j3 - Ji -j2,m 1 +m 2 + m 3 \J - j 2 ,m 2 + m 3 ,ji,m 1 )(j 3 - j 2 ,m 2 + m 3 \j 2 ,m 2 , J,m 3 ) = 

(44) 

7.3 Wigner D-Matrix 

The irreducible representation of SO (3) are called Wigner D-matrices T> e g and are indiced by an 
integer £ — 0, . . . , oo. The Ah order representation works on a C 2£+1 dimensional vector space. We 
denote the components of T) e g by D^^g). In Euler angles in ZYZ-convention we have 

Z?L(7,/3,«) = e- im ^L(/3)e inQ , (45) 
where d e mn (/3) is the 'small' Wigner d-matrix, which is real-valued and explicitly written as 

d e mn (p) = [(^ + m)!(£-m)!(^ + n)!(£-n)!] 1 / 2 V, 7 - — ( ,7 1)m "" + * n 

mn\^ J LV 1 /V / V / \ / J L^s {£+n—s)\s\(m—n+s)\(£—m—s)\ (Af\\ 

x (cos f ) (sin f ) 

The representations of different order are connected via the Clebsch Gordan coefficients by: 

DL= Dt^D^ilmlhmiJ^ilnlhn^hn^ (47) 

mi +mo=m 



m^+m2=m 
n^-\-n2 —ti 



and 

D^nAn 2 = E ^L(^Ki^i^2m 2 )(/n|/im,/ 2 n 2 ) (48) 



Another important equality is 



/ dg D[ lk Di, m ,Di m = ^^(jnlfn'Jk^ijmlfm'Jk) (49) 



7.4 Mixed Quadratic Terms 

The action of the terms T T J±i and J±iT T can be computed directly from the equations ([7]), ([8]), 
([11]) to 

(T T J±if)L = "i E E VW + l)/2-(m±l)(m±2)/2^^ 

j /= j-1,j,j+1 n=n'-\-q 

g=-l,0,l 

(HjV, lg)0m|/(m ± 1), 1( T 1)) 0^ m±2 
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and 



(J±iT T f) J nm = Vj(j + l)/2-m(m±l)/2 

j ,= j — + 1 n=n' + q 



g=-l,0,l m =m / ^2 



7.5 Proof of equation (14) 



From equation ([32]) we directly know the corresponding equality for the spherical derivatvies: 

O 7 _L_ i 

J,M 1 

Setting f i = 4 = 1 we know that the sum over J takes only values for J = and J = 1. By 
rotating the frame of reference d e ~D e (g) T d e we get 

2J + 1 

t 



(Di(^) T 9 1 ) mi (D 1 (^) T 9 1 ) m2 = — ^Hlmxllro* JM)<10|10, JOXD'fo)"^)* 



J,M 



with 7^ = (D 1 (g) d ) m and evaulting the Clebsch Gordan coefficients we end up with 



10 2 



r mi r m2 = T - V E (i^i|i^2,2M) (d 2 (s) t £> 2 ) m 



M=-2 

which was to show. 
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